ENERGY STABLE AND MOMENTUM CONSERVING HYBRID 
FINITE ELEMENT METHOD FOR THE INCOMPRESSIBLE 
NAVIER STOKES EQUATIONS 
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Abstract. A hybrid method for the incompressible Navier-Stokes equations is presented. The 
method inherits the attractive stabilizing mechanism of upwinded discontinuous Galerkin methods 
when momentum advection becomes significant, equal-order interpolations can be used for the ve- 
locity and pressure fields, and mass can be conserved locally. Using continuous Lagrange multiplier 
spaces to enforce flux continuity across cell facets, the number of global degrees of freedom is the same 
as for a continuous Galerkin method on the same mesh. Different from our earlier investigations on 
the approach for the Navier-Stokes equations, the pressure field in this work is discontinuous across 
cell boundaries. It is shown that this leads to very good local mass conservation and, for an appro- 
priate choice of finite element spaces, momentum conservation. Also, a new form of the momentum 
transport terms for the method is constructed such that global energy stability is guaranteed, even 
in the absence of a point-wise solenoidal velocity field. Mass conservation, momentum conservation 
and global energy stability are proved for the time-continuous case, and for a fully discrete scheme. 
The presented analysis results are supported by a range of numerical simulations. 
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1. Introduction. A method that combines attractive features of discontinuous 
and continuous Galerkin hnite element methods for the incompressible Navier Stokes 
equations was presented in Labeur and Wells [l[ and further extended to the case of 
moving domains and free-surface flows in Labeur and Wells Q ■ The method incorpo- 
rated naturally the evaluation of upwinded advective fluxes on cell facets, in the same 
spirit as discontinuous Galerkin methods, thereby stabilizing flows with significant 
momentum advection, and it is possible to use equal-order polynomial bases for the 
velocity and pressure components. However, the number of global degrees of freedom 
on a given mesh is the same as for a continuous Galerkin method using the same poly- 
nomial orders. The issue regarding the significantly greater number of global degrees 
of freedom for low- to moderate-order discontinuous Galerkin methods compared to 
continuous Galerkin methods is thus circumvented. However, the method in [l[ was 
restricted to continuous pressure fields, and it could not be proven that the method 
is globally energy stable. These short-comings are addressed in this paper, with a 
formulation presented that permits discontinuous pressure fields, is globally energy 
stable, conserves momentum and has excellent local mass conservation properties. 

The key to the methodology that we present for constructing schemes is the 
postulation of cell-wise balances, subject to weakly enforced boundary conditions. 
The boundary condition to be satisfied (weakly) is provided by a function that lives 
on cell facets only. An equation for this extra field is furnished by insisting on weak 
continuity of the associated 'numerical' flux. The concept of weak enforcement of 
flux continuity across cell facets is central in hybridized finite element methods (for 
an overview see Q). A feature of these methods is that functions on cells are linked 
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to functions on neighboring cells only via functions that live on cell facets, and not 
directly via the flux terms. Therefore, functions on cells can be eliminated locally 
in favor of the functions that live on cell facets only (via static condensation), thus 
reducing the number of globally coupled degrees of freedom. If the functions enforcing 
the continuity of the fluxes, and which live only on cell facts, are discontinuous, then 
point-wise continuous fluxes can be obtained for suitably chosen function spaces. In 
contrast, in our method we advocate the use of facet functions that are continuous, 
which leads to a significant reduction in the number of globally coupled degrees of 
freedom, since the local elimination procedure will lead to a global problem of the 
same size as a corresponding continuous Galerkin method. Yet, it is straightforward to 
demonstrate local momentum and mass conservation, in terms of the numerical fluxes, 
as is typical of discontinuous Galerkin methods. Also the stabilizing mechanism of 
the flux formulation, involving the advection terms and the pressure- velocity coupling, 
are directly inherited from the discontinuous Galerkin method and lead to favorable 
stability properties. 

We are not alone in considering methods that draw on both discontinuous and 
continuous Galerkin methods. Hughes et al. J4J developed a method for the advection- 
diffusion equation, and the formulation in [l| for the advection-diffusion equation re- 
duces to that of Hughes et al. Q in the advective limit. In the diffusive case, there 
is a subtle difference, with the diffusive flux in Hughes et al. Q being upwinded, 
whereas a centered approach is used in Labeur and Wells Simulations using the 
approach for the advection-diffusion equation exhibited very good stability proper- 
ties, minimal dissipation and standard convergence rates. For the case of the linear 
advection-diffusion-reaction equation, stability (via an inf-sup condition) and con- 
vergence at a rate of k + 1 in the diffusive limit and k + 1/2 in the advective limit 
was later proved In the context of hybridized methods, Cockburn and co-workers 
have published a number of works (e.g. @, Q) that share features with the method- 
ology that we consider. A hybrid field on cell interfaces is presented in Egger and 
Schoberl Q for the advection-diffusion problem. Giizey et al. [§] present a hybrid 
continuous-discontinuous finite element method for elliptic problems, coined embed- 
ded discontinuous Galerkin method, which is conceptually related to the method in 
Labeur and Wells [l[ . 

The method formulated and analyzed in this work is an extension of the method 
presented in Labeur and Wells [l| for the advection-diffusion equation and for the 
incompressible Navier-Stokes equations. Unlike in our previous efforts [H, Q, we con- 
sider here pressure fields that are discontinuous across cell facets. The impact of this 
on the local (cell-wise) mass conservation properties of the method is demonstrated. 
Another feature that distinguishes the formulation developed in this work from our 
earlier work for the Navier-Stokes equations is the use of a skew-symmetric form of 
the advective term. The derivation of the skew-symmetric formulation is not trivial 
in the considered setting, but it is shown that it brings the advantage of guaranteeing 
stability in terms of the total kinetic energy, even when the velocity field is not point- 
wise solenoidal. The combination of discontinuous pressure fields and skew-symmetric 
advection terms leads to a method that for equal-order basis functions preserves mass 
and momentum and is also stable in terms of the total kinetic energy. The analy- 
sis results that we present are supported by a number of computations for both the 
Stokes and incompressible Navier-Stokes equations. The computer code necessary to 
reproduce all examples presented in this work is available in the supporting material 
[lOj under a GNU public license. 
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The remainder of this work is structured as follows. We first define concretely 
the problem of interest, and then develop a semi-discrete finite clement formulation. 
Some properties of the semi-discrete problem are then analyzed. This is followed by a 
particular fully-discrete formulation, and it is shown that the considered properties of 
the semi-discrete problem are inherited by the fully discrete problem. This is followed 
by numerical examples, after which conclusions are drawn. 

2. Incompressible Navier Stokes equations. We consider a domain of in- 
terest fl C R d , where 1 < d < 3 is the spatial dimension. The boundary dfl is 
assumed sufficiently smooth and the outward unit normal vector on dfl is denoted by 
n. The boundary is partitioned such that L^nLjv = dfl and L/jULat = 0. The time 
interval of interest is I — (0, tjv]. 

The non-dimensional incompressible Navier-Stokes problem on fl x / reads: given 
the viscosity v, the forcing term / : !lx/-> R d , the momentum flux h : IV x / — > M. d 
and the solenoidal initial condition u^: fl — >• R d , find the velocity field u: fix I — > R d 
and the pressure field p : fl x I — > R such that 

OH 

— + V-<r = / on fix I, (2.1) 

er =pl -2vV s u + u®u on fix I, (2.2) 

V • u = on fl x I, (2.3) 

u = on T D x /, (2.4) 

ern — max (u ■ n, 0) u = h on Tji x /, (2-5) 

u (x, 0) = uo (x) on fl, (2.6) 

where a is the momentum flux, I is the identity tensor, \7 s u = (Vit + Vit T ) /2 is 
the symmetric gradient, in which [Vm]^ = dui/dxj, and [u ® u]ij = UiUj . The 
Neumann boundary condition has been formulated such that on portions of IV on 
which u ■ n < (inflow boundaries) the total momentum flux is prescribed, while on 
portions of IV on which u ■ n > (outflow boundaries) only the diffusive part of the 
momentum flux is prescribed. 

3. Finite element method. The hybrid finite element method is defined in 
this section. The essence of the method is posing all balance equations cell-wise in a 
weak sense, with a suitably constructed numerical flux, and complementing the cell- 
wise balance laws by a global equation enforcing weak continuity of the numerical flux 
across cell facets. 

3.1. Definitions. We consider a triangulation T of the domain fl into open, non- 
overlapping sub-domains K (cells). The outward unit normal vector on the boundary 
dK of each cell is denoted by n. Adjacent cells share a common facet F, and T — (J F 
is the union of all facets, including the exterior boundary facets. A measure of the 
size of a cell K is denoted by \%k- When evaluated on a shared facet, Kk is used to 
imply the average cell size measure of the adjacent cells. 

Consider first the vector finite element spaces Vh and Vh- 

V h := {v h e [L 2 (T)] d ,v h G [P k {K)] d Vifer}, (3.1) 

V h := [v h e [L 2 (F)] d , v h G [P~ k (F)] d V F G T, v h = on T D } , (3.2) 

where Pk(K) denotes the space of Lagrange polynomials on K of order k > 0, and 
Pk{F) denotes the space of Lagrange polynomials on F of order k > 0. The space Vh 
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contains vector-valued functions that are discontinuous across cell boundaries, while 
functions in Vh are defined on cell facets only. Furthermore, functions in Vh satisfy 
the homogeneous Dirichlet boundary condition on Tjj. Scalar finite element spaces 
Qh and Qh are defined by: 

Q h ■= {q h e L 2 (T) , q h e P m (K)VKeT}, (3.3) 

Qh ■= {Qh G L 2 (JT) , q h G P rn (f)VfeJ}, (3.4) 

where the polynomial orders m > and fh > 0. Mirroring the vector spaces, Qh 
contains functions that are discontinuous across cell facets, while functions in Qh are 
defined on cell facets only. 

For algorithmic reasons, it may be advantageous to compute with the finite ele- 
ment spaces 

V^:=V h n[H 1 (J=)] d , (3.5) 
Qt :=Q h nH\F), (3.6) 

in place of Vh and Qh, respectively, using polynomial orders k > 1 and fh > 1. This 
will be discussed in Section [6j and all computational results presented in Section [7] 
will employ facet functions that are continuous. 

3.2. Semi-discrete weak local/global balances. We formulate now a semi- 
discrete finite clement problem by considering what we will refer to as 'local' and 
'global' equations. The 'local' equations solve the problem cell-wise in which the 
velocity and pressure boundary conditions are provided by auxiliary fields that live 
on cell facets only. To determine the fields that live on cell facets, 'global' equations 
are formulated by requiring weak continuity of the mass and momentum fluxes across 
element interfaces. The methodology behind the construction of the formulation is 
elucidated by presenting a collection of Galerkin problems for the various balances, 
after which the considered Galerkin finite element problem is completely and formally 
defined. 

3.2.1. Local/global continuity equation. A Galerkin approximation of the 
incompressibility constraint (|2.3p in a cell- wise fashion requires that the approximate 
velocity it/,. £ Vh satisfies 

^ / u h ■ Vq h dx-y^ l Uh -n q h ds = V q h E Q h , (3.7) 

K JK k JdK 

where iih is the 'numerical' mass flux on dK, and is chosen to be 

u h = u h - —^r (ph ~ Ph) n, (3.8) 
v + 1 

in which ph £ Qh and ph G Qh are pressure fields and j3 > is a parameter required 
for stability when using equal-order basis functions for the velocity components and 
pressure fields. When using a lower polynomial order for the pressure field relative to 
the velocity field it is possible to use j3 = [ll[. Penalization of the pressure jump was 
used by Hughes and Franca [l2| to stabilize equal-order methods with discontinuous 
pressure for the Stokes equation, and by other authors for discontinuous Galerkin 
methods IH 14 1. However, different from Hughes and Franca [T^, we add a non- 



dimensional unit viscosity to the term in the denominator to permit consideration 
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of the inviscid limit. Central in p.8[) is that the pressure-stabilizing term involves 
the difference between ph and ph, rather than the jump in ph across a facet as in 



other works 12], LL3J, [l4[ . Equation (I3.7P is local in the sense that there is no direct 
interaction between ph on neighboring cells. This is a key feature of the method with 
practical implications that will be elaborated upon in Section [5] 

The numerical mass flux in (|3.8p is not unique on cell facets; it can take on different 
values on different sides of a facet. This is in contrast with standard discontinuous 
Galerkin methods, in which the numerical mass flux is constructed such that it is 
uniquely defined on facets. A 'global' continuity equation is now furnished by insisting 
that the numerical mass flux iih be weakly continuous across cell facets, in that it 
satisfies 



E 

K 



/ ii h ■ n q h ds - I u h ■ n q h ds = V qh £ Qh, (3.9) 
JdK Jan 



where Uh £ VJ,. Note that (|3.9p implies that iih ■ n = Uh ■ n (weakly) on d£l. 

3.2.2. Local/global momentum balance in conservative form. At time t 
and given the forcing term / £ [L 2 (T)] d , the viscosity v, the velocity iih £ Vh, and 
pressures ph £ Qh and ph £ Qh, consider a Galerkin approximation Uh £ Vh that 
satisfies the following weak formulation of the momentum balance (|2.1[) : 



/ ~~57~ ' v h dx - J~] [ <r h : Vv h dx + f cr h n ■ v h ds 

Jn dt ^J K ^ JdK 

+ V / 2v(u h -Uh)-V s v h nds = f v h dx Vv h £V h , (3.10) 
K JdK Jn 

where the momentum flux <Th on cells is given by 

cr h = a {u h ,Ph) = Phi - 2vV s u h +u h ® u h , (3.11) 

and the 'numerical' momentum flux &h on cell boundaries is given by 

= &a,H + &&,Hi (3.12) 

where the advective flux & a ,h is 

o- a ,h = &a (u h ,u h ,Ph,Ph) = u h ®ii h + (u h - u h ) ® \u hl (3.13) 



in which iih is given by (|3.8p . A is a function that takes on a value of either one or 
zero and is defined below, and the diffusive flux &d,h is 

crd.h = crd (Uh, u h ,Ph) = Phi - 2v\/ s u h - t—^v i u h - u h ) ® n. (3.14) 

tlx 

The function A takes on a value of one on inflow cell boundaries (where u n < 0) , and 
takes on a value of zero on outflow cell boundaries (where u n > 0). The formulation 
for the advective interface flux involves upwinding since <r a = u ® u on outflow cell 
boundaries and rr a = u ® u on inflow cell boundaries. In (|3.10p . the fourth term on 
the left-hand side ensures symmetry of the diffusion operator (see Arnold et al. [15]). 
In (|3.14p . a is a penalty parameter, and such a term is typical of interior penalty 
methods. Just as for standard interior penalty methods, the role of the penalty term 
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in this context is to ensure stability, as detailed in Wells Equation (|3.10[) is 'local' 
in the sense that the weak momentum balance equation is posed cell-wise. 

As with the numerical mass flux Q3.8P , the numerical momentum flux Q3.12[) is not 
single- valued on cell facets. A 'global' momentum balance equation is therefore formu- 
lated by insisting on continuity of the numerical flux er/,. across element facets. This 
continuity constraint is imposed weakly by requiring that, for a given flux boundary 
condition h 6 [L 2 (T N )] d , 

y2 &hn-Vhds= / (1 — A) (uh ® H>h) n • ds 
K JdK Jr N 

+ [ h-v h ds Vv h e V h . (3.15) 
Jr N 

The above equation implies that the numerical momentum flux (T^n and the momen- 
tum flux on r^r, given by h + (1 — A) (uh ® Uh) n, coincide in a weak sense. 

3.2.3. Local/global momentum balance in advective form. We now re- 
phrase the conservative forms of the local and global momentum balance equations 
into advective formats with a view to formulating a skew-symmetric version of the 
advective terms. 

Considering first the local momentum equation, substitution of the fluxes (13. . 
(|XT2j) and (j3"T3l into (jBTTUl) yields 

/ ~KT -Vhdx-y [ a d ,h ■ Vuft da; - V" / {u h ® u h ) : Vv h dx 
Jn ot k-* k it 

+ T2 Vd,hn -v h ds + y2 {u h - n) u h ■ v h ds 

K JdK K JdK 

+ ^2 \{ii h ■ n) (u h - u h ) ■ v h ds 
K JdK 

+ / 2v(uh-u h )-'V s Vhnds= / f-v h dx, (3.16) 
^ JdK Jn 



in which <Td,h = Phi — 2v'SI s ui l is the diffusive flux on cells. Applying partial integra- 
tion to the advective terms on K in (|3.16p . 

/ - v hdx-y~] [ a- d ,h ■ Vu/j dx + y^ [ {Vu h u h ) ■ v h dx 
Jn at % Jk k J k 

+ T2 ( v ' u h) u h ■ v h dx + y~] / ((u h - u h ) ■ n) u h ■ v h ds 

K J K K JdK 

+ y~] I CTd^n -v h ds + y2 x ( u h ■ n) (u h - u h ) ■ v h ds 

K JdK K JdK 

+ V / 2v{u h -Uh)-V s v h nds= / f-v h dx. (3.17) 
„ JdK Jn 



We choose to discard the integrals involving V • Uh and (v,h ~ **/,) -n, which by virtue 
of the continuity equation (|3.7[) and under an appropriate regularity assumption on 
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the exact solution will not disturb consistency of a Galerkin scheme (this will be 
addressed formally in Section [4]). A reduced version of (|3 . 1 T[> now reads: 

/ -Vhdx + y] [ (Vtth u h ) ■ v h dx - V* / a d . h : Vv h dx 
Jn ot k-* k it 

+ V" / \(u h ■ n) (u h - u h ) ■ v h ds + & dih n-v h ds 

K JdK K JdK 

+ V / 2v(u h -u h )-V s v h nds= / f-v h dx. (3.18) 
K JdK Jn 

Considering next the global momentum balance, inserting the expressions for the 
numerical flux (I3.12p and (|3.13[) into the global flux continuity equation (|3.15p yields, 



y~] / trd,h.n -Vh ds + y] / (u h ■ n) u h ■ v h ds 

K JdK K JdK 

+ A (v,h ■ n) (u h - uh) ■ v h ds 

K JdK 

- (l-X)(u h -n)uh-v h ds= / h-v h ds. (3.19) 
The second integral in (I3.19[) can be expanded as 

V] / {u h -n)u h - v h ds = V] / (u h ■ n) (u h - u h ) ■ v h ds 

K JdK K JdK 

+ V" / {{uh - Uh) ■ n) u h ■ Vh ds + / (u h ■ n)u h ■ v h ds, (3.20) 
K JdK Jdn 

where we have used that is single-valued on cell facets, by definition. Discarding 
the term involving (v,h — Uh) ■ n, which is consistent with continuity equation (|3.9p . 
and substituting (I3.20[) into (|3.19[) leads to the following advective form of the global 
momentum equation: 

V] / &d,hn -v h ds-y2 (1 - A) (u h ■ n) (u h - u h ) ■ v h ds 

K JdK K JdK 

j A (u h ■ n) u h ■ v h ds = j h- v h ds, (3.21) 



where we have taken into account that Dj, = on <9S!\r 



JV- 



3.3. Semi-discrete finite element formulation. We define now a collection 
of functionals that together will define a complete finite element problem. For conve- 
nience, the notation U '.= (u,u,p,p) will be used. 

Based on the local continuity equation (|3.7[) . we define the functional: 



F C (U; q) :— 2_] / u-\7qdx — 2_, / u nqds, 
K Jk k JdK 



(3.22) 
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where F c is linear in q. Similarly, from the global continuity equation (|3.9p . the 
functional 



F C (U; q) ■= 2_. / unqds— I u-nqds, 
JdK Jon 



(3.23) 



is defined, where F c is linear in q. For the momentum equations, we define a local 
momentum balance functional that is a weighted combination of the local conservative 
balance (|3.16p and the local advective balance (|3.18p : 



F m (U; v) :— f ^ • v dx — \ / ( u ® u) : \7v dx 
Jn ot — J K 

+ (1 — x) 5^ / u ) ' v + X / (u ■ n) u ■ vds 

K Jk k JdK 

+ y^ / A (it ■ n) (u — u) ■ v ds — y^ / er^ : Vi> da; 

K Jok K Jk 

+ V / & d n-vds + y2 2v(u-u) -V s vnds- I f-vdx, (3.24) 

K JdK K JdK Jo. 

where x € [0; 1] an d F m is linear in v. In the same fashion, the global momentum 
flux continuity equations (|3.19p and (|3.21l) are weighted and summed, leading to, 



F m (U; v) :— x ^ / (u ■ n) u ■ v ds — (1 — x) /J / (it ■ n) (u — u) ■ v ds 

K JdK K JdK 

+ 2_. / X(u • n) (u — u) ■ vds + yj / &d.n • * ds 

— / (x — A) (m • n) m • w ds — / fa, •■yds, (3.25) 

where .F m is linear in v. 
Defining now 

F(U; W) := F m (U; v) + F m {U; v) + F C {U; q) + F C (U; q), (3.26) 

where W — (v,v,q,q), a semi-discrete finite element problem at time t involves: 
given the forcing term / £ [L 2 (f2)] d the boundary condition fa G [L 2 (r^)] d and the 
viscosity v, find U \ S Vh X Vh x Qh x Qh such that 

F(E/ h ; W fc ) = V W h £V h xV h x Qh x Qh- (3.27) 

This completes the formulation of the semi-discrete finite element problem. 

4. Properties of the semi-discrete formulation. We now demonstrate the 
consistency, mass conservation, momentum conservation and energy stability proper- 
ties of the method for the semi-discrete formulation in (|3.27p . The presented results 
hold for the spaces Vh and Qh, and deliberately also for the more restrictive case 
in which Vh and Qh are replaced by V h * and Q* t , respectively, which we advocate in 
practice and will use in numerical examples. 

Proposition 4.1 (consistency). If at a given time t, u € (H 2 (Q)\ and p e 
H 1 (fi) solve equations (|2.1l) - (|2.5p . and u = j(u) and p = j(p) on T . where 7 is a 
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trace operator, then 

F(U; W h ) = V W h EV h xV h x Q h x Q h , (4.1) 

for any \ S [0, 1]. 

Proof. Considering first Vh = 0, Vh = and qh — 0, applying integration by parts 
to (|4~T|) leads to 



V f (V • u) q h dx - V / — ^- (p - p) q h ds = 



= Vq h eQ h , (4.2) 



which holds due to satisfaction of (|2.3p and p = 7(p). Setting t)/j = 0, Vh — and 
q h = in fn}, 

V] / u-nq h ds+ u-nq h ds- u-nq h ds = V q h & Qh, (4.3) 
^ JdK\dn Jan Jan 



which holds due to the regularity of u and because u = ^{u). 

Setting Vh = 0, qh = and q~h = in (|4.ip and applying integration by parts, 

(jjj; + V • o- - /J • v h dx - (1 - x) jf (V • u) w • dx 

+ x / ((** — u) • ti) tt • U/j ds + 7^ / X (u • n) (u — u) • Vh ds 

K JdK K JdK 

+ / (P ~ P) 71 ' Vh ^ s — y2 / -^—2v{u-u)®n-v h ds 

K JdK K JdK K 

+ V / 2^(it-«)-V s v, l nds = Vi)^^, (4.4) 

which holds due to u and p satisfying equations (|2.1j) and (|2.3I) . the regularity of u 
and because p = j(p) and w = u = 7(11). Finally, setting Vh = 0, qh = and qh =0 
in 

/ 

IdK K JdK u JdK 

ds 



/ ern-Vhds + 'S^ / (p — p) n • Vh ds — S~] / ——2v{u — u)®n-Vhds 

K JdK K JdK K JdK "-K 

+ / ((it — u) ■ n) u ■ Vh ds + / X (u ■ n) (u — u) ■ Vh< 

K JdK K JdK 

- (1 - x) I {u ■ n)u ■ v h ds - (x- X)(u-n)u- v h ds 

K JdK Jt n 

= [ h-v h ds VveVh, (4.5) 
Jr N 

which holds due to the regularity of w, because p = 7(p) and ft = u = j(u) and due 
to satisfaction of the flux boundary condition (|2.5I) . Equation (|4.1I) follows from the 
summation of (I4.2I) - (|4.5I) and the linearity of F in v, v, q and q. □ 

Key to the proof of Proposition 14. II is the consistent formulation of the numerical 
fluxes, that is <r = cr and ii = u if u = u and p = p. 

Proposition 4.2 (mass conservation) . If Uh and Uh satisfy (|3.27p . then 

u h -nds = VKeT (4.6) 

dK 
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u h -nds = 0. (4.7) 

Proof. Setting Vh — Vh = 0, qh = 1, and qt = 1 on the cell K and qh = on 
T\K leads to (gU). Setting v h = v h = and q h = q h = 1 in (j3~27)) leads to (pLT]) . □ 
The local conservation property is in terms of the numerical mass flux it, as is 



typical for discontinuous Galerkin methods applied to Stokes flow [13j, uM, ■ Classical 
local mass conservation would be satisfied if (3 — 0, but (3 must be greater than zero 
for stability of the equal-order case [l2j]. If the pressure field is chosen to be one 
polynomial order lower than the velocity field, then (3 can be set equal to zero and 
mass is conserved locally and exactly. However, it will be shown that reducing the 
size of the pressure space requires compromising on either momentum conservation 
or energy stability. 

Proposition 4.3 (momentum conservation) . If Uh and Uh solve (13.271) . and the 

function spaces Vh, Vh, Qh and Qh are selected such that for a constant but otherwise 
arbitrary vector c it holds that Vh ■ c G QhS Vh £ Vh and Vh ■ c e QhS Vh £ Vh, then 

<l_ 

di 



4 / u h dx = [ fdx - [ & h nds V KeT, (4.8) 

' J K JK JdK 



and ifTr> = 



^- [ u h dx= [fdx- f (1 - A) (uh ■ n) u h ds - f hds. (4.9) 
"* Jn Jn Jen Jon 



Proof. Setting Vh — ej and qh = — (1 — x) u h • &j on K, where ej is a canonical 
unit basis vector, Vh = and qh = on T\K, Vh = and q~h = in (|3.27[) . 

^- { u h ■ e dx + [ (u h ■ n) u h ■ e 3 ds + [ A (u h ■ n) (u h - u h ) ■ ej ds 

dt JK JdK JdK 

+ / cr d . h n ■ e 3 ds = / f-ejdx, (4.10) 

JdK J K 

which from the definition of the fluxes in equation (|3.13[) proves (|4.8[t . 

Sett ing v h = ej, v h = -ej, q h = -(1 - X) u h • e 3 and q h = -(I - x)uh • e 3 
in ([3~27| leads to (]S9j directly. □ 

Local momentum conservation is in terms the numerical flux &h, as is a typical 
feature of discontinuous Galerkin methods. Note also the requirement on the size of 
the pressure space relative to the components of the velocity space, which would not 
be satisfied by methods that use lower-order polynomials for the pressure than for 
the velocity, such as Taylor-Hood elements. Such elements only conserve momentum 
when conservative forms of the momentum equation are used which, in the advective 
limit, requires compromising on energy stability. Provided that the requirements on 
the sizes of the function spaces are met, momentum conservation holds irrespective 
of the value of x, i.e. for advective as well as conservative forms of the advection 
operator, and it will be shown that for \ = 1/2 the method is also energy stable (see 
proposition I4.4[) . 

For cases with Dirichlet boundary conditions (Tjj ^ 0), demonstration of momen- 
tum conservation is less straightforward since Vh can not be set equal to ej on Tp- 
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This difficulty can be overcome by introducing auxiliary flux terms on T]j. Details of 
the approach can be found in Ref . [l6| ■ 

PROPOSITION 4.4 (global energy stability). If u h solves (|3.27p with \ = 1/2 
and homogeneous boundary conditions, then in the absence of forcing terms and for 
a suitably large a 



d_ 

dt 



\u h \ 2 dx < 0. (4.11) 



Proof. Setting v h — u h , v h — -Uh, qh = ~Ph and Qh = ~Ph m P-271) gives, for 
X = 1/2, 

J ■ u h + J Q - A^j (uh ■ n) \u h - u h \ 2 ds 

- / (?d,h ■ Vm/, dx - V" / a- d ,h n ■ - Uh) ds 

K JK K JdK 

+ I 2v(u h -u h )-V s u h nds+ [ ( ^ - A ) (u h ■ n) \u h \ 2 ds 
K JdK Jr N \ z / 

- I u h ■ Vpu dx - V] / Uh ■ n (p h ~ Ph) ds + u h -np h ds = 0. (4.12) 

K J K K JdK JdQ. 



Substituting the expressions for the diffusive fluxes given in (|3.11[) and (|3. 14|) into the 
preceding equation, 

/ ^ ~5T dx + X / ^ I"' 1 ' n l \ Uh - u h \ 2 ds + ( 2v |V s m^| 2 dx 
Jn z Ct J gK Z k ^ K 

+ / T~~ Zv \ u h — u h\ 2 ds + 2 / 2v {uh — Uh) ■ V s un ds 

K JdK h K JdK 

- V" / Phi- V« fi dx - V] I p h n ■ (u h - u h ) ds + I - \u h ■ n\ \u h \ 2 ds 

K Jk k JdK Jr N Z 

- V" / u h ■ Vph dx - / u h ■ n(p h - p h ) ds + u h -np h ds = 0, (4.13) 

K Jk k JdK JdQ. 



in which we have also used (1/2 — A) u ■ n = \u ■ n\ /2 on facets. After substitution 
of the mass flux Uh given in Q3.8P and the application of partial integration of the 
pressure gradient term, we finally obtain, 

4: \ \ \ u hf dx + V" / \ \u h ■ n\ \u h - u h \ 2 ds + V] f 2v \V s u h \ 2 dx 
at Jo. 1 x JdK z x Jk 

+ T/ j—Zv \uh - u h \ 2 ds + / 2v (u h - u h ) ■ V s u h nds 

~ JdK h K JdK 

+ f l -\u h -n\\u h \ 2 ds + V / ^L\p h - Ph f ds = 0. (4.14) 
Jr N z K JdK i + v 

For the case v — 0, all terms in (I4.14j) other than the time derivative term, are 
guaranteed to be non-negative, and therefore (|4.11l) holds. For the case v > 0, no 
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conclusion as to the sign of (uh — Uh) ■V s ui l n in (|4.14l) can be drawn. However, there 
exists an a > 0, independent of Yik, such that (|4.1ip holds (see the proof in [B| for 
the diffusion equation). □ 

The key to the energy stability is the use of the combined conservative/advective 
form of the momentum equations. The total kinetic energy in the method will de- 
crease monotonically, despite Uh not being point-wise divergence-free. The amount 
of dissipation is determined by the difference between the cell and facet velocity fields 
on facets (iih — Uh) and the difference between the cell and facet pressure fields on 
facets (ph — Ph)- It is not possible to prove a cell- wise kinetic energy inequality. 

Recall that the energy stable scheme is also momentum conserving if the di- 
mensional components of the velocity space are subspaces of the pressure space (see 
Proposition I4.3|) . Otherwise, simultaneous momentum conservation f (|4.8p and (|4.9|1 ) 
and energy stability (I4.11[) is not possible. Notably, this will be the case for finite 
element methods using lower-order polynomials for the pressure than for the veloc- 
ity. For advection dominated flows such elements can only guarantee energy stability 
when compromising on momentum conservation or by adding artificial viscosity. In 
the viscous limit the requirement on the size of the velocity and pressure spaces can be 
relaxed without compromising energy stability, permitting a lower polynomial order 
for the pressure space than for the velocity space, which for Stokes flow is advanta- 
geous from the viewpoint of accuracy, as will be shown in Section 17.11 

A complete stability proof, for the Stokes case alone, would require a more subtle 
analysis, with the usual stability conditions demonstrated for suitably defined norms 
that include both functions on cells and the functions on facets. A priori stability 
and convergence estimates for the method applied to the advection-reaction-diffusion 
equation have been proved [5], and efforts in this direction for the Stokes equations 
are ongoing. 



5. A fully-discrete formulation. We now present a fully-discrete formulation. 
The time interval of interest, /, is partitioned such that I = (0, t%, . . . , £at_i, and 
time increments are denoted 5t n — t n+ \ — t n . We consider a ^-method for dealing 
with the time derivative, with mid-point values of a function y given by 



y n+e ■= (1 - 9)y n + 0y n +i, (5.1) 



where 9 G [0, 1] is a parameter. 

The advective velocity will be evaluated at the current time t ni thereby lineariz- 
ing the problem (Picard linearization). For the momentum-related F-functionals pre- 
sented in Section 13.31 we now present time-discrete counterparts. The term A is 
always evaluated on the basis of the known velocity field at time t n . A time-discrete 
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counterpart of p. 241) reads 

Fst,m(U n+ i;v) := / — ^di-xV / (u n+g ® « n ) : Vt> da; 

+ (1 - / (Vun+eM,,) • u da; + \ V] / (tt„ • n) u n+e • w ds 

^ ^ ./9K 

+ X(u n -n) (u n +e - u n+g ) ■ vds - V] / <x d .„ +( > : Vi> da; 

+ a-d.n+en ■ vds + 2v (u n+e - u n+e ) -V s vnds 

K JdK K JdK 



/ fn+e - vdx i ( 5 - 2 ) 



and a time-discrete counterpart of (|3.25p reads 
Fst,m(Un+i;v) ■= X^ / (m„ • n) tx n+e • w ds 

- (1 - x) V" / n) (ttn+e - • * ds 

~ J dK 



+ ^2 X(ii n ■ n) (u n+g - u n+e ) ■ vds + V] / & d:n+g n ■ v ds 

K JdK K JdK 

- (X - X) (u n ■ n) u n+e ■ v ds - \ h n+e -vds. (5.3) 

JFn J^n 

Defining 

Fst(U n+1 ; W) := F St . m {U n+1 ;v) 

+ F st!m (U n+1 ;v) + F c (U n+1 ;q) + F c (U n+1 ; q), (5.4) 

a fully-discrete finite element problem at time t n+ \ involves: given the solution Uh.n £ 

Vh x Vh x Qh x Qh at time t n , the forcing term f n+g G [L 2 (fi)] d , the boundary 

condition h n+ g 6 [L 2 (1^)]^ and the viscosity v, find Uh, n +i £ V/, X t4 X Qh X 
such that 

F 5i (t/ h)n+1 ; = V W h e V* x V), x Q h x (5.5) 

For the scheme that has been adopted, F$ t is linear in both Uh, n +i and Wh- 

We now demonstrate that the considered fully discrete formulation inherits the 
conservation and energy stability properties of the semi-discrete case. As for the semi- 
discrete case, all results hold if the spaces Vh and Qh are replaced by V h * and Q^, 
respectively. 

Proposition 5.1 (fully discrete mass conservation). If Uh, n +\ an d Uh, n +i sat- 
isfy (|5.5|) . then 

u h , n+ i ■ n ds = V K e T, (5.6) 

dK 
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u h ,n+i ■ nds = 0. (5.7) 

an 



Proof. The proof follows the same steps as the proof of Proposition 14.21 □ 
Proposition 5.2 (fully discrete momentum conservation). IfUh,n and Uh, n +i 
solve (|5.5I) . and the function spaces Vh, Vh, Qh and Qh are selected such that for 
a constant but otherwise arbitrary vector c it holds that Vh ■ c G QhS Vh G Vh and 
Vh- c G QhS v h G Vh, then 

f Uh >"+\- Uh > n dx = [ f n+e dx- [ & h ,n + ends VKeT, (5.8) 

JK ot JK JdK 

and ifTr) = % 

Uh,n+1 — Uh,- 



St ^-JJ^d, 

- (1 - A) (u h , n ■ n) u h „n+e ds- h n+g ds. (5.9) 
Jan Jan 

Proof. The proof follows the same steps as the proof to Proposition 14.31 Equa- 
tion (I5.8[) follows from setting = ej and qh = — (1 — x) u h.n+e • ej on K, Vh = 
and qh = on T\K, Vh = and q~h = in (|5.5p . Equation (15.91) follows from setting 
v h = ej, v h = -e^ q h = -(1 - x)uh,n+e ■ e 3 and q h = -(1 - x)uh,n+e ■ ej in ([53]). □ 

Proposition 5.3 (fully discrete energy stability). If Uh,n and Uh,n+i solve 
(|5.5[) with x — 1/2, f — and Tu = dQ, for 9 > 1/2 and suitably large a 

/ \u h ,n+i\ 2 dx < / |xt^,„| 2 dx. (5.10) 

Proo/. For x = 1/2, setting v h = u h , n+ e, v h = -Uh,n+e qh = -9p h , n +6 and q h = 
-0ph,n+e in (l53|) . and adding to ([53)1 F c (U n ; -(l-9)p n+e ) and F c (U n ; -(l-9)p n+ e), 

J ttfe -"+ 1 ^ . it h n+9 dx + '^2 J Q - A^ (Mh,n • n) |«h, n +e - w/j.n+e | 2 ds 

I <7d,h,n+e ■ Vmi,„ + « dx - >J / o-d.h,n+en ■ (u h , n +e - M/i,n+e) ds 

+ ^2 ^ V (™h.n+9 ~ Uh.n+e) ' Vllft^+Ofl ds 

x ./air 

+ y Q - A^j • n) lii/j^i+el 2 ds - ^ J u h , n +e ■ ^Ph.n+e dx 

I Uh,n+e ■ n(p h ,n+8 - Ph.n+e) ds + u h ,n+e-nph, n +eds = 0. (5.11) 
Taking into account that 

u h<n+e =St[9--\ + - (5.12) 
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leads to 



Uh.n+l - U h n ( 1\ f \u h n+ i - U h , n \ 2 

u h , n+e dx = - - / ■ da; 



< > st V - / ./<> st 

2 1 r U.. |2 



l /■ K n+1 | d i /■ l^ da; _ (513) 



2 in 5i ' 2 J n St 

Following the same steps as in Proposition ^. 41 proves that (I5.10[) holds when 9 > 1/2 
and for sufficiently large a. □ 



6. Algorithmic aspects. The fully-discrete finite element problem in f|5.5[) can 

lead to an efficient numerical implementation in which the functions on cells (v,h t n+i 
and ph^n+i) are eliminated cell- wise in favor of the functions that live only on cell 
facets {uh, n +i and ph, n +i) via static condensation. Key to this algorithmic feature 
is that functions on cells are not linked directly across cell facets, in contrast with 
conventional discontinuous Galerkin methods. Rather, functions on neighboring cells 
communicate via the functions defined only on cell facets. Moreover, if the functions 
Uh and ph are chosen to be continuous, the method will result in the same number 
of global degrees of freedom as for a continuous method on the same mesh (if interior 
degrees of freedom are eliminated from a continuous method via static condensation) . 
Despite having the same number of global degrees offreedom as a continuous method, 
stabilizing mechanisms that are typical of discontinuous Galerkin methods are natu- 
rally incorporated. Moreover, for the advection-diffusion equation, it has been proved 
that the approach has the same stability properties as classical upwinded discontin- 
uous Galerkin methods [f|. More detailed discussions on algorithmic aspects can be 
found in Refs. 0,0. 

7. Examples. We present now a number of examples in support of the analysis 
presented in the preceding sections. The computer code used to compute the exam- 
ples presented in this section has been generated automatically from expressive input 



using tools from the FEniCS Project (http : //www . f enicsproject . org) 17]. Specif- 
ically, an expressive domain-specific language for finite element variational statements 
in combination with automated code generation IH 19|, 2(| and a programmable en- 



vironment has been used [17J. The computer code used for all examples presented in 
this work is available under a GNU public license in the supporting material [l(| ■ 

All examples use triangular elements, uniform partitionings and continuous facet 
functions, that is Uh 6 V£ and ph 6 Q* h - When computing errors, analytical solutions 
which are polynomial are represented exactly. Otherwise analytical solutions are 
interpolated using eighth-order Lagrange basis functions on the same mesh. Exact 
quadrature is used in all cases. 

7.1. Stokes flow with source. We consider a Stokes problem (by neglecting 
the momentum advection terms) with v = 1 on a unit square with u = on dtt. The 
source term / is chosen such that the exact solution is: 



u x - x 2 (1 - x y (2y - 6y 2 + 4y 3 ) , 

u y = —y 2 (1 - yf (2.x - 6x 2 + Ax 3 ) , (7.1) 
p = x (1 — x) . 

The constraint J^pdx = 1/6 is enforced by means of a Lagrange-multiplier, matching 
the solution in (17. 1[) . 



16 



R. J. LABEUR AND G. N. WELLS 



We investigate convergence rates in the L 2 -norm for the pressure and the velocity 
fields using equal-order polynomial elements (k = k = m = in). Polynomial orders 
ranging from one to five are considered. These results complement those presented in 
Labeur and Wells [1] for the same boundary-value problem, but in which the pressure 
field was continuous and only a polynomial order of one was considered. We set 
a = 6k 2 , based on observations for higher-order elements Q, and use /3 — 10 -4 . The 
observed convergence behavior is presented in Figure I7TT1 Standard convergence rates 
of order k + 1 for the velocity field and of order k for the pressure field arc observed. 
The error in the divergence of the velocity field is examined via 



and the computed ediv is shown in Figure 17.21 for various polynomial orders and h- 
rcfincmcnt. Clearly, the divergence error is small. For comparison, the divergence 
errors using the same method for the velocity field, but with a continuous pressure 
field [l| are shown in Figure 17.31 The observed convergence rates are similar to 
those for the discontinuous pressure case. However, particularly for the lower-order 
elements, the divergence error is significantly greater in the continuous pressure case. 

For comparison, we consider a Taylor-Hood element with a continuous piecewise- 
quadratic velocity field and a continuous piecewise-linear pressure field, and an ele- 
ment constituted of a continuous piecewise-quadratic velocity field, enriched by cubic 
bubble functions, and a discontinuous piecewise-linear pressure field. The latter ap- 
proach is referred to by some authors as the Crouzeix-Raviart method (e.g. [2l|, HH), 
a convention that we adopt here. For our method we use corresponding polynomial 
orders of k = k = 2 for the velocity and m = m = 1 for the pressure and penalty 
parameters a — 6k 2 and (3 = 0. Recall that it is permitted to use f3 = in this case 
since the polynomial degree of the pressure field is lower than the polynomial degree 
of the velocity field, see Section l3.2.1l The observed convergence behavior is presented 
in Figure 17.41 which shows the expected convergence rates for the Taylor-Hood and 
Crouzeix-Raviart methods, and with the method formulated in this work showing the 
same rates, which are also the same as for the k = k = m = m = 2 case presented in 
Figure EU The divergence error, measured by ediv, is shown in Figure 1731 We note 
that while the Crouzeix-Raviart method conserves mass locally, the divergence error 
when measured in ediv is very close to that of the Taylor-Hood method, whereas for 
our method the divergence error is effectively zero. 

These convergence results demonstrate that in the viscous limit the equal-order 
pressure approximation is sub-optimal from the viewpoint of accuracy. However, in 
advection dominated flows the simultaneous satisfaction of Proposition 14.31 (momen- 
tum conservation) and Proposition 14.41 (energy stability) relies on a sufficiently rich 
pressure field relative to the velocity field, i.e. k < m, k < fh, and consequently j3 > 0. 

7.2. Kovasznay flow. We now consider the incompressible Navier-Stokes equa- 
tions by examining the following analytical solution due to Kovasznay j23[: 




(7.2) 



■it 



'X 



■y 



1 - e Xx cos (27ry) , 
^e Xx sin {2ny), 



(7.3) 
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where C is an arbitrary constant and 



_ Re 



Re 1 



4tt' 



1/2 



(7.4) 



where Re is the Reynolds number. The solution represents laminar flow in the wake 
of a grid, see also 24, 25} . 

We use a rectangular domain 51 := {(x,y) £ (—0.5,1) x (—0.5,1.5)}. On 9f2 





Dirichlet boundary conditions for the velocity are specified according to equation (|7.3[) . 
The pressure is prescribed in the lower-left corner of the domain. Equal-order poly- 
nomial elements are used (k = k = m = m) with polynomial orders ranging from 
one to five. The parameters \ — 1/2, a — 6k 2 and j3 = 10 -4 are used. We solve the 
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(a) 



10 



10 



a, 10 -3 
I 

a, 



10 



10 




io- 



10 



10° 



(b) 

Fig. 7.4. Stokes flow: computed L 2 errors in (a) velocity and (b) pressure with h-refinement 
for polynomial orders k = k = 2 and m = m = 1 (a = 6fc 2 and /3 = 0), compared with Taylor-Hood 
(T-H) and Crouzeix-Raviart (C-R) methods. 



stationary problem using a fixed point iteration with stopping criterion 



el +1 + e* 



< TOL, 



(7.5) 



where e l u and e„ +1 are the L 2 velocity error norms, relative to the exact solution, of 
the consecutive iterates i and i + 1, respectively, and TOL is a given tolerance which 
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10- 2 




h 



Fig. 7.5. Stokes flow: divergence error with h-refinement for polynomial orders k = k = 2 and 
m = m = 1 (a = 6fc 2 and /3 = 0), compared with Taylor-Hood (T-H) and Crouzeix-Raviart (C-R) 
methods. 



is set to 10~ 4 . 

For Re — 40, the observed convergence rates in the L 2 -norm for the velocity 
and pressure fields are presented in Figure 17.61 Convergence rates of order k + 1 for 
the velocity field and of order k for the pressure field are observed. It was verified 
that these convergence results also hold for x — (advective scheme) and x = 1 
(conservative scheme). 

7.3. Backward-facing step flow. The next example concerns stationary two- 
dimensional flow over a backward-facing step. Figure 17.71 presents the set-up of the 
problem. The step height S is equal to half the height of the main channel height D 
and the velocity profile in the inflow channel is parabolic with maximum velocity J7 max . 
Behind the step a recirculation zone develops with the re-attachment length xt de- 
pending on the Reynolds-number, as investigated experimentally by Armaly et al. 
[26j j . The Reynolds number is defined as 

Re=^-, (7.6) 

where U is two-thirds of the maximum inflow velocity and v is the kinematic viscos- 
ity 0. 

The numerical test concerns the comparison of experimental and computed values 
of the dimensionless reattachment length xi/S [26| . We consider zero step length (L = 
0), and a rectangular computational domain f2 := {(x,y) 6 (0,15) x (0,1)} which 
extends from the step over a length of 30 times the step height in the downstream 
direction. The domain is partitioned using 301 vertices in ^-direction and 31 vertices in 
y-direction. On the left boundary (x = 0) the parabolic velocity profile with C/ max = 1 
is imposed for 1/2 < y < 1 using a Dirichlct boundary condition. Along the outflow 
boundary (x — 15) a homogeneous Neumann boundary condition for the velocity is 
used. On all other boundaries u = 0. The pressure degree of freedom for p in the 
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Fig. 7.8. Backward- facing step: comparison of measured (•) and computed reattachment lengths 
for polynomial orders of k = 1 (solid) and k = 2 (dashed), experimental data from \2&1 . 



lower-left corner of the domain is fixed. Using the definition in (|7.6I) . the maximum 
inflow velocity and kinematic viscosity are adjusted to obtain a range of Reynolds 
numbers between 100 and 800. Equal order polynomials are used (k = k = m — fa) 
with k — 1 and k = 2 considered. The parameters \ — 1/2, a — 6k 2 and f3 = 10 -4 are 
used. We solve the stationary problem using a fixed point iteration with a stopping 
criterion based on the L 2 -norm of the velocity, 



Uh 



|4+1 

lo,n 



\Uh 



Q.Q 



\U h 



lo,n 



< TOL, 



(7.7) 



where i and i + 1 denote successive iterates and TOL is a tolerance, which is set 
to 10~ 6 . 

The computed dimensionless reattachment lengths for various Reynolds num- 
bers are presented in Figure 17.81 and are compared against measured data [26j . For 
Re < 400 the computed results are in good agreement with the results obtained from 
the experiments. For Re > 400 the computed results gradually deviate from the 
measurements, in a similar way as the results computed by Kim and Moin [27j ]. and 
which can be attributed to the emergence of three-dimensional flow structures [26| . 
The computed streamlines for Re — 800 and polynomial orders of one are shown in 
Figure 17.91 The computed streamlines involve a secondary recirculation bubble which 
resides between dimensionless distances of 10.4 and 20.1 from the step , which is in 
good agreement with observed values of 11.2 and 19.6, respectively 26 1. 



7.4. Chaotic advection. We now consider the energy stability properties of the 
method for the incompressible Navier-Stokes equations with x — 1/2- The incom- 
pressible Navier-Stokes equations are solved on the unit square with zero viscosity 
and boundary conditions u ■ n = and h ■ s = on 9f2, where s ■ n = 0. This 
corresponds to impermeable free-slip boundaries. The pressure is prescribed to be 
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Fig. 7.10. Relative change in total kinetic energy between time steps for the incompressible 
Naviei — Stokes test with \ = 1/2, 9 = 1/2 and v = for linear (k = 1 ) and quadratic (k = 2) 
elements. In both cases, a mesh with 32 X 32 vertices is used. 



zero at a point in the domain. The initial condition tto = is used. A time step 
5t = 0.2 is adopted and the mesh has 32 cell vertices along each axis. To create a 
chaotic velocity field, in the first simulation step a random forcing term / is used. 
Uniform random variables are generated at vertices such that for each component of 
the forcing vector fa 6 [—1,1]. This field is then interpolated using linear Lagrange 
finite element basis functions. For the first step, v = 1 x 1CP 5 , after which it is set 
to zero. This is done to start the simulation, since with Uq = 0, if v = then (|5.5|) 
cannot be solved as the advective velocity at t = is zero. For the first 5 time steps, a 
backward Euler scheme is used (9 = 1) to damp oscillations due to the discontinuous 
nature of the forcing term. After the first 5 steps, 9 = 1/2 is used. 

The relative change in the total kinetic energy between steps once 9 — 1/2 is 
presented in Figure [7.101 for the case of linear basis functions (fe = 1) for all fields 
and the case of quadratic basis functions (k = 2) for all fields. Consistent with the 
analysis, the kinetic energy is observed to decrease monotonically. Not unexpectedly, 
the relative dissipation is smaller for the k = 2 case. 

8. Conclusions. A generalization of a hybrid method that inherits attractive 
properties of continuous and discontinuous Galerkin methods has been presented and 
analyzed for the incompressible Navier-Stokes equations. The method incorporates 
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upwinding of the advective momentum flux naturally, it is observed to be stable for 
equal-order velocity/pressure basis functions and it has very good local mass conser- 
vation properties. These properties, usually associated with discontinuous Galerkin 
methods, can be achieved with the same number of global degrees of freedom as a 
continuous Galerkin method on the same mesh, thereby obviating the common crit- 
icism of discontinuous Galerkin methods that the number of degrees of freedom is 
too large compared to continuous methods. In contrast with our earlier work, the 
new formulation presented here involves a pressure field that is discontinuous across 
cell facets. This has implications for local mass conservation, which in the presented 
formulation is guaranteed in terms of the numerical flux. It is shown that with ap- 
propriately chosen (equal order) function spaces the method conserves momentum. 
Moreover, the new formulation presented in this work uses a skew-symmetric form of 
the momentum advection term. It has been shown that this, in combination with a 
suitable time integration scheme, guarantees that the global kinetic energy will de- 
cay monotonically, even if the velocity field is not point-wise divergence-free. The 
properties of the method that have been demonstrated by analysis are supported by 
numerical examples. Standard convergence rates for a range of polynomial orders are 
observed in the Stokes and Navier-Stokes examples, and simulations comparing the 
continuous and discontinuous pressure cases illustrate the advantage of discontinuous 
pressure fields for local mass conservation. The Navier-Stokes example concerning 
the flow over a backward facing step shows that the method performs well in an 
advection dominated case. The Navier-Stokes example concerning the evolution of 
a randomly generated velocity field demonstrates the energy decaying property of 
the skew-symmetric momentum advection term. The complete computer code for 
performing all presented numerical examples is made freely available under an open 
source license as part of the supporting material. 
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